
/*
 ============================================================================
 Name        : validateMatches.c
 Author      : Petko Fiziev
 Version     :
 Copyright   : Your copyright notice
 Description : Validate matches from the blind search
 ============================================================================
 */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <ctype.h>


#include "DIGEST_V2/SuffixTree.h"
#include "DIGEST_V2/ExactSearch.h"

#include "utils.h"

extern int maxOccurences;

typedef struct {
	long id;
	int real_length;	
	WORD *bitseq;
	int count;
	ResultRecord *occurences;	
} SEQUENCEMATCH;

//static int compare_sequences(const SEQUENCE *seq1, const SEQUENCE *seq2) {
static int sort_by_file(const void *pseq1, const void *pseq2) {
	SEQUENCEMATCH *seq1 = *(SEQUENCEMATCH * const *)pseq1;
	SEQUENCEMATCH *seq2 = *(SEQUENCEMATCH * const *)pseq2;
	
	if (seq1->count == 0) return -1;
	if (seq2->count == 0) return 1;
	
	int fdiv = seq1->occurences[0].fileNumber - seq2->occurences[0].fileNumber;
	if (!fdiv) {
		return seq1->occurences[0].startPos - seq2->occurences[0].startPos;
	} else {
		return fdiv;
	}
}

static int sort_by_id(const void *pseq1, const void *pseq2) {
	//SEQUENCEMATCH *seq1 = *(SEQUENCEMATCH * const *)pseq1;
	//SEQUENCEMATCH *seq2 = *(SEQUENCEMATCH * const *)pseq2;
	int id1 = (*(SEQUENCEMATCH * const *)pseq1)->id;
	int id2 = (*(SEQUENCEMATCH * const *)pseq2)->id;
	if (id1 > id2) { return 1; }
	else if (id1 < id2) { return -1; } 
	else { return 0; }
		
}




int main(int argc, char *argv[]) {

	int i;
	char seqfilename[MAX_PATH_LENGTH] = {0};

		
	if(argc<6) {
		fprintf(stderr, "To run: %s <treefolder> <treefileprefix> <blindSearch_output> <original_fasta_file> <maximum_occurences>\n", argv[0]);
		return 1;
	}
	
	maxOccurences = atoi(argv[5]);
	int max_mismatches = 0;
	char buf[MAXLINE] = {0};
	char charseq[MAXLINE] = {0};
	
	FILE *input = fopen(argv[3], "rb");
	if (input == NULL) {
		fprintf(stderr,"cannot open %s\n", argv[3]);
		return 1;
		}
	
	// read total sequences, total bitseq length and total occurences
	int total_sequences, total_count, total_bitseq_length;
	fseek(input, -3 * sizeof(int), SEEK_END);
	int intbuf[3];
	if (fread(intbuf, sizeof(int),3, input) != 3) {
		fprintf(stderr,"cannot read seq info from %s\n", argv[3]);
		return 1;
		}
	total_sequences = intbuf[0];
	total_bitseq_length = intbuf[1];
	total_count = intbuf[2];
	rewind(input);
	
	fprintf(stderr,"total sequences: %d\n", total_sequences);	
	fprintf(stderr,"total bit seq length: %d\n", total_bitseq_length);	
	fprintf(stderr,"total occurences: %d\n", total_count);	
	
	SEQUENCEMATCH **seqs = (SEQUENCEMATCH **)malloc(total_sequences * sizeof(SEQUENCEMATCH *));
	SEQUENCEMATCH *seqdata = (SEQUENCEMATCH *)malloc(total_sequences * sizeof(SEQUENCEMATCH));
	WORD *bitseqdata = (WORD *)calloc(total_bitseq_length, sizeof(WORD));
	ResultRecord *occurences = (ResultRecord *)malloc(total_count * sizeof(ResultRecord));
	
	
	if (seqs == NULL || seqdata == NULL || bitseqdata == NULL || occurences == NULL){
		fprintf(stderr, "Not enough memory!\n");
		}
	
	
	fprintf(stderr, "Reading input | ");
	start_watch();	
		
	
	int bitseqlen = 0;
	int bitseqdata_index = 0;
	int occurences_index = 0;
	for (i = 0; i < total_sequences; i++) {
		seqs[i] = &seqdata[i];
		
		if(fread(&seqs[i]->id, sizeof(long), 1, input) != 1) {
				fprintf(stderr,"cannot read seq id from %s\n", argv[3]);
				return 1;
			}				
		// read the real sequence length
		if(fread(&seqs[i]->real_length, sizeof(int), 1, input) != 1) {
			fprintf(stderr,"cannot read real_length from %s\n", argv[3]);
			return 1;
		}
		// read the bit sequence
		bitseqlen = blen(seqs[i]->real_length);
		seqs[i]->bitseq = &bitseqdata[bitseqdata_index];
		bitseqdata_index += bitseqlen;
		if(fread(seqs[i]->bitseq, sizeof(WORD), bitseqlen, input) != bitseqlen) {
			fprintf(stderr,"cannot read from %s\n", argv[3]);
			return 1;
		}
		// read the number of occurences
		if(fread(&seqs[i]->count, sizeof(int), 1, input) != 1) {
			fprintf(stderr,"cannot read count from %s\n", argv[3]);
			return 1;
		}
		seqs[i]->occurences = &occurences[occurences_index];
		
		int to_read = MIN(seqs[i]->count, maxOccurences);
		occurences_index += to_read;
		// read the file numbers
		if (to_read > 0) {
			if(fread(seqs[i]->occurences, sizeof(ResultRecord), to_read, input) != to_read) {
				fprintf(stderr,"cannot read occurences from %s\n", argv[3]);
				return 1;
			} 
		}
			
		//bitseq2char(seqs[i]->bitseq, seqs[i]->real_length, buf);
		//printf(">%ld\tcount:%d\n%s\n", seqs[i]->id, seqs[i]->count, buf);
					
	
	}
	stop_watch();	
	
	fprintf(stderr, "Sorting by match positions| ");
	start_watch();	
	
	qsort(seqs, total_sequences, sizeof(SEQUENCEMATCH *), sort_by_file);
	stop_watch();	
	
	
	char outputfilename[MAXLINE];
	sprintf(outputfilename, "%s.validateMatches", argv[3]);
	FILE *output = fopen(outputfilename,"wb");
	if(output == 0) {
		fprintf(stderr,"can't open file: %s.\n", outputfilename);	
		return 1;
	}	

	
	start_watch();	
	fprintf(stderr, "Validating matches | \n");
	fprintf(stderr, "reading sequence files: ");
	FILE *seqfile = NULL;
	int current_seq_file = -1;
	int plen = -1;
	int k = -1;
	int mismatches = 0;
	for (i = 0; i < total_sequences; i++) {		
		if (seqs[i]->count > 0) {
			
			if (current_seq_file != seqs[i]->occurences[0].fileNumber) {
				if (seqfile != NULL) { fclose(seqfile); }
				current_seq_file = seqs[i]->occurences[0].fileNumber;
				//fprintf(stderr, "csf: %d\n", current_seq_file);
				sprintf(seqfilename,"%s%s_%d", argv[1], argv[2], current_seq_file);
				fprintf(stderr, "%d, ", current_seq_file);
				if ((seqfile = fopen(seqfilename,"r")) == NULL) {
					fprintf(stderr, "Cannot open sequence file: %s\n", seqfilename);
					return 1;
					}
				}
			plen = seqs[i]->real_length;
			//fprintf(stderr, "plen:%d\n", plen);
			
			fseek(seqfile, seqs[i]->occurences[0].startPos, SEEK_SET);
			//fprintf(stderr, "fseek: %d\n", seqs[i]->occurences[0].startPos);
			if ( fread(buf, sizeof(char), plen, seqfile) != plen) {
				fprintf(stderr, "Cannot read from: %s\n", seqfilename);
				return 1;
			}
			//fprintf(stderr, "fread: %s\n", buf);
			
			bitseq2char(seqs[i]->bitseq, seqs[i]->real_length, charseq);
			/*
			fprintf(stderr, "\ns: %s\n", charseq);
			fprintf(stderr, "b: %s\n", buf);
			fprintf(stderr, "count: %d\n", seqs[i]->count);
			for (k = 0; k < seqs[i]->count; k++) {
				fprintf(stderr, "\tANGS_fn[%d]:%d\tANGS_startPos[%d]:%d\n", i, 
															seqs[i]->occurences[k].fileNumber, 
															i, 
															seqs[i]->occurences[k].startPos);
			} */
			mismatches = 0;
			for(k = 0; k < plen; k++) {
				//if(tolower(charseq[k]) != tolower(buf[k])) {
				if(charseq[k] != buf[k]) {
					mismatches++;
					if (mismatches > max_mismatches) {
						seqs[i]->count = 0;
						}
				}
			}	
		}
		
		//bitseq2char(seqs[i]->bitseq, seqs[i]->real_length, buf);
		//printf(">%ld\tcount:%d\n%s\n", seqs[i]->id, seqs[i]->count, buf);
		
		/*
		fprintf(output,">%ld\tcount:%d", seqs[i]->id, seqs[i]->count);
		for (k = 0; k < seqs[i]->count; k++) {
			fprintf(output, "\tfn:%d\tstartPos:%d", seqs[i]->occurences[k].fileNumber, seqs[i]->occurences[k].startPos);
		}
		fprintf(output, "\n%s\n", charseq); */
	}	
	fprintf(stderr," done | ");
	stop_watch();	
	
	
	fprintf(stderr, "Sorting by id | ");
	start_watch();	
	qsort(seqs, total_sequences, sizeof(SEQUENCEMATCH *), sort_by_id);	
	stop_watch();	
	
	
	fprintf(stderr, "Writing output | ");
	start_watch();	
	FILE *fasta_input = fopen(argv[4], "r");
	if (fasta_input == NULL) {
		fprintf(stderr,"cannot read from %s\n", argv[4]);
		return 1;
		}
		
	for (i = 0; i < total_sequences; i++) {
		fseek(fasta_input, seqs[i]->id, SEEK_SET);
		if(!fgets(buf, MAXLINE, fasta_input)) {
			fprintf(stderr,"cannot read from %s\n", argv[4]);
			return 1;
		}
		plen = strlen(buf);
		if (buf[plen-1] == '\n') { buf[plen-1] = 0; }
		
		bitseq2char(seqs[i]->bitseq, seqs[i]->real_length, charseq);
		int to_output = MIN(maxOccurences, seqs[i]->count);
		fprintf(output,"%s\tANGS_total:%d\tANGS_output:%d", buf, seqs[i]->count, to_output);
		
		for (k = 0; k < to_output; k++) {
			fprintf(output, "\tANGS_fn[%d]:%d\tANGS_startPos[%d]:%d", k, 
															seqs[i]->occurences[k].fileNumber, 
															k, 
															seqs[i]->occurences[k].startPos);
		}
		fprintf(output, "\n%s\n", charseq);
	}
		
	
	stop_watch();	
	
	fclose(fasta_input);
	fclose(input);
	fclose(output);
	return 0;
}
